%% Codes for the article "Dissection Mechanisms of Financial Crises: Intermediation and Sentiment", by Arvind Krishnamurthy and Wenhao Li.
%% Set Path and Add Path (not needed if we call this file from "Main_Step1.m")
% clear all;
% mainpath = fileparts(mfilename('fullpath'));
% addpath(genpath([ mainpath, '/Codes_Model' ] ))
% cd(mainpath);

S1 = load( 'solutions/baseline_benchmark.mat' );
S2 = load( 'solutions/baseline_rational.mat' ) ; 
S3 = load( 'solutions/baseline_behavioral.mat' );


%% ------------------- Load Diagnostic Solution ---------------------------------------------
model_sol = post_processing(S3.model_sol);
par = S3.model_sol.par;
w_dim = kron(ones(par.N_lambda,1), par.w_grid );
lambda_dim = kron(par.lambda_grid, ones(par.N_w,1));
xK_matrix_fun =  scatteredInterpolant(w_dim,   lambda_dim ,  reshape( model_sol.xK_matrix, par.N_w*par.N_lambda, 1 ),   'linear')  ;
bank_credit_GDP_fun = scatteredInterpolant(w_dim,    lambda_dim ,  reshape( model_sol.psi_matrix .* model_sol.p_matrix  ./ model_sol.productivity_matrix , par.N_w*par.N_lambda, 1 ),   'linear')  ;
TB_loaded = readtable( 'simulated data/lambda_relationship_diagnostic.xlsx');
lambda_diagnostic = TB_loaded.lambda_diagnostic;
lambda_rational = TB_loaded.lambda_rational;
index_w = 35;
w = par.w_grid(index_w);
xK_rational = xK_matrix_fun( w*ones(size(lambda_rational)),  lambda_rational );
xK_diag = xK_matrix_fun( w*ones(size(lambda_rational)),  lambda_diagnostic );
bank_credit_rational = bank_credit_GDP_fun(w*ones(size(lambda_rational)),  lambda_rational );
bank_credit_diag = bank_credit_GDP_fun(w*ones(size(lambda_rational)),  lambda_diagnostic );


%% ------------------------- Figure 3b: Diagnostic lambda v.s. rational lambda, average relationship from simulation -------------------------
disp('--------------------------------------------------------------------------------------------------------')
disp('--------------------- Figure 3b: Diagnostic Lambda v.s. Rational Lambda -------------------------------------------------------')
figure( 'Position', [0 0 330 320] ) ;  
% fit the average relationship
slm_fitted = slmengine( lambda_rational,  lambda_diagnostic);
lambda_diag_fitted = slmeval(lambda_rational, slm_fitted);
plot( lambda_rational, lambda_rational,  lambda_rational, lambda_diag_fitted , '--', 'LineWidth', 1  ); xlim([0,0.7]);
lgd = legend( 'Bayesian belief', 'Diagnostic belief' );  % xlabel('Rational lambda (expected intensity of liquidity shocks)'); ylabel('Belief of liquidity shock intensity');
xlabel('$\lambda^{Bayesian}$', 'Interpreter','latex');  ylabel('$\lambda^{belief}$', 'Interpreter','latex'); 
legend boxoff; lgd.FontSize = 10; lgd.Location='northwest'; 
set(gca,'FontSize',lgd.FontSize);
saveTightFigure( [figure_path ,'/model_illustration/diagnostic_against_rational_lambda.pdf']  );



%% ---------------------- Figure 5: Credit spread against bank loan spread and liquidity premium --------------------
disp('----------------------------------------- Figure 5 ---------------------------------------------------------------')
% Figure 5(a)
disp('--------------------- Figure 5(a): Liquidity Premium and Credit Spread -------------------------------------------------------')
model_sol = post_processing(S2.model_sol, true);
par = model_sol.par;
w_index = [40, 45];
rf_rd_diff =  model_sol.lambda_matrix .* par.alpha ./ ( 1- model_sol.xK_matrix.*model_sol.kappap_matrix  - par.alpha*model_sol.delta_x_matrix   );
figure( 'Position', [0 0 360 320] ) ;  
plot( model_sol.credit_spread(w_index(1),:),  rf_rd_diff(w_index(1),:), model_sol.credit_spread(w_index(2),:),  rf_rd_diff(w_index(2),:), '--' , 'LineWidth', 1.5  ); 
xlabel( 'credit spread'  ); ylabel( 'liquidity premium'  );   legend( {['w=',num2str(round(par.w_grid(w_index(1)),2))],['w=',num2str(round(par.w_grid(w_index(2)),2))]}, 'Location', 'northwest' )
saveTightFigure( [figure_path ,'/model_illustration/credit_spread_against_liq_prem.pdf']  );
% Figure 5(b)
disp('--------------------- Figure 5(b): Loan Spread and Credit Spread -------------------------------------------------------')
excess_returns =   (model_sol.sigmap_matrix + par.sigmaK).^2.*model_sol.xK_matrix +  model_sol.lambda_matrix .* model_sol.kappap_matrix ./ ( 1- model_sol.xK_matrix.*model_sol.kappap_matrix  - par.alpha*model_sol.delta_x_matrix   )  - model_sol.lambda_matrix .* model_sol.kappap_matrix;
excess_returns_rd = excess_returns + rf_rd_diff;
figure( 'Position', [0 0 360 320] ) ;  
slm1=slmengine( model_sol.credit_spread(w_index(1),:), excess_returns_rd(w_index(1),:) , 'concavedown','on' );
slm2 = slmengine(model_sol.credit_spread(w_index(2),:),  excess_returns_rd(w_index(2),:), 'concavedown', 'on');
plot( model_sol.credit_spread(w_index(1),:), slmeval(model_sol.credit_spread(w_index(1),:), slm1) , model_sol.credit_spread(w_index(2),:),  slmeval(model_sol.credit_spread(w_index(2),:), slm2) , '--' , 'LineWidth', 1.5  ); 
xlabel( 'credit spread'  ); ylabel( 'bank loan spread'  );   legend( {['w=',num2str(round(par.w_grid(w_index(1)),2))],['w=',num2str(round(par.w_grid(w_index(2)),2))]}, 'Location', 'northwest' )
saveTightFigure( [figure_path ,'/model_illustration/credit_spread_against_capital_excess_return.pdf']  );


%% ---------------------- Figure 10: Belief and bank credit, average relationship from simulation ----------------------------
disp('--------------------------------------------------------------------------------------------------------')
disp('------------------------- Figure 10: Belief and Bank Credit ------------------------------')
figure( 'Position', [0 0 360 320] ) ;  
slm_fitted = slmengine( lambda_rational,  bank_credit_diag);
bank_credit_diag_fitted = slmeval(lambda_rational, slm_fitted);
plot( lambda_rational, bank_credit_rational, 'b',  lambda_rational,  bank_credit_diag_fitted, 'r--', 'LineWidth', 1  ); xlim([0,0.7]);
lgd = legend( 'Bayesian belief', 'Diagnostic belief' ); xlabel('Rational lambda (expected intensity of liquidity shocks)'); ylabel('bank credit/GDP');
legend boxoff; lgd.FontSize = 10;
set(gca,'FontSize',lgd.FontSize);
saveTightFigure( [figure_path , '/model_illustration/belief_and_bank_credit.pdf']  );


%% ---------- Figure 16: bank credit and equity return in the benchmark model ------------------------
disp('--------------------------------------------------------------------------------------------------------')
disp('------------------------ Figure 16: Bank Credit and Equity Return in the Benchmark Model ------------------------------')
model_sol = post_processing(S1.model_sol, true);
figure( 'Position', [0 0 400 300] ) ;  
xval = mean(model_sol.xK_matrix .* model_sol.w_matrix ./ model_sol.productivity_matrix,2);
yval = mean(  - model_sol.kappaw_matrix + model_sol.muw_matrix*1/12,  2);  
plot(xval, yval, 'LineWidth', 1.5);
xlabel('bank credit/GDP');  ylabel('one-month % change of bank equity after dN shock')
saveTightFigure( [figure_path, '/model_illustration/bank_credit_and_equity_return_Benchmark.pdf']  );